#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>

#include <glew.h>
#include <GL/glut.h>

#include "tetrahedral.h"
#include "gpudt/pba3D.h"


#define iWidth 16
#define iHeight 16
#define iDepth 16

//*************** input demo********************//
void randomPoints(int texSize, int nPoints);
void demo1(int texSize, int nPoints); //2*2*2 + 8
void demo2(int texSize, int nPoints); //2*2*2 + 7
void cubeDemo(int texSize, int nPoints);
//**********************************************//

//**********************************************************************************************
//* Input parameters
//**********************************************************************************************
/*int fboSize     = 64;
int nVertices   = 100;*/ 

int fboSize     = 256;
//generateRandomPoints
//int nVertices   = 50;
//int nVertices   = 2*2*2 + 8; //demo 1
//int nVertices   = 2*2*2 + 7; //demo 2
//int nVertices = 300;
//int nVertices   = 2*2*2 + 2;
int nVertices   = 300763;
void (*demoFun)(int texSize, int nPoint) = cubeDemo;
int faceColor = 1;
int display = 1;

int phase1Band  = 1; 
int phase2Band	= 1; 
int phase3Band	= 2; 

#define TOID(x, y, z, w)    ((z) * (w) * (w) + (y) * (w) + (x))

typedef struct {
    double totalDistError, maxDistError; 
    int errorCount; 
} ErrorStatistics; 

//**********************************************************************************************
//* Global Vars
//**********************************************************************************************
int *inputPoints, *inputVoronoi, *outputVoronoi; 
Tetrahedral *tetrahedral = NULL;
BoundaryMissing *boundaryMissing = NULL;
int boundaryCount = 6;
int tetraCount = 0;
ErrorStatistics pba; 

GLuint texture; 
float scale = 0.01f, rotateX = 0.0, rotateY = 0.0; 
bool isLeftMouseActive = false, isRightMouseActive = false; 
int oldMouseX = 0, oldMouseY = 0; 

//**********************************************************************************************
//* Random Point Generator
//**********************************************************************************************
// Random number generator, obtained from http://oldmill.uchicago.edu/~wilder/Code/random/
unsigned long z, w, jsr, jcong; // Seeds
void randinit(unsigned long x_) 
{ z =x_; w = x_; jsr = x_; jcong = x_; }
unsigned long znew() 
{ return (z = 36969 * (z & 0xfffful) + (z >> 16)); }
unsigned long wnew() 
{ return (w = 18000 * (w & 0xfffful) + (w >> 16)); }
unsigned long MWC()  
{ return ((znew() << 16) + wnew()); }
unsigned long SHR3()
{ jsr ^= (jsr << 17); jsr ^= (jsr >> 13); return (jsr ^= (jsr << 5)); }
unsigned long CONG() 
{ return (jcong = 69069 * jcong + 1234567); }
unsigned long rand_int()         // [0,2^32-1]
{ return ((MWC() ^ CONG()) + SHR3()); }
double random()     // [0,1)
{ return ((double) rand_int() / (double(ULONG_MAX)+1)); }

void randomPoints(int texSize, int nPoints) {
	int tx, ty, tz, id; 
	//randinit((unsigned long)time(NULL));
	randinit(0);
	for (int i = 0; i < nPoints; i++)
	{
		do { 
			tx = int(random() * texSize); 
			ty = int(random() * texSize); 
			tz = int(random() * texSize); 
			id = TOID(tx, ty, tz, texSize); 
		} while (inputVoronoi[id] != MARKER); 

		inputVoronoi[id] = ENCODE(tx, ty, tz); 
		inputPoints[i] = ENCODE(tx, ty, tz);  
		//printf("\n%d, %d, %d", tx, ty, tz);
	}
}

void demo1(int texSize, int nPoints) {
	//int nVertices   = 2*2*2 + 8;
	int tx, ty, tz, id; 
	int nCount = 0;
	int offset = 10; 
	int step = 40;
	int size = 64;
	int len = (size - 2*offset) / step + 1;
	for(int i=0; i<len; i++) {
		tx = i*step + offset;
		for(int j=0; j<len; j++) {
			ty = j*step + offset;
			for(int k=0; k<len; k++) {
				tz = k*step + offset;
				id = TOID(tx, ty, tz, texSize); 
				inputVoronoi[id] = ENCODE(tx, ty, tz); 
				inputPoints[nCount++] = ENCODE(tx, ty, tz);
			}
		}
	}

	tx = 32; ty = 32; tz = 63;
	id = TOID(tx, ty, tz, texSize); 
	inputVoronoi[id] = ENCODE(tx, ty, tz); 
	inputPoints[nCount++] = ENCODE(tx, ty, tz);

	tx = 0; ty = 63; tz = 63;
	id = TOID(tx, ty, tz, texSize); 
	inputVoronoi[id] = ENCODE(tx, ty, tz); 
	inputPoints[nCount++] = ENCODE(tx, ty, tz);

	tx = 63; ty = 63; tz = 63;
	id = TOID(tx, ty, tz, texSize); 
	inputVoronoi[id] = ENCODE(tx, ty, tz); 
	inputPoints[nCount++] = ENCODE(tx, ty, tz);

	tx = 0; ty = 32; tz = 32;
	id = TOID(tx, ty, tz, texSize); 
	inputVoronoi[id] = ENCODE(tx, ty, tz); 
	inputPoints[nCount++] = ENCODE(tx, ty, tz);

	tx = 32; ty = 63; tz = 32;
	id = TOID(tx, ty, tz, texSize); 
	inputVoronoi[id] = ENCODE(tx, ty, tz); 
	inputPoints[nCount++] = ENCODE(tx, ty, tz);

	tx = 63; ty = 32; tz = 32;
	id = TOID(tx, ty, tz, texSize); 
	inputVoronoi[id] = ENCODE(tx, ty, tz); 
	inputPoints[nCount++] = ENCODE(tx, ty, tz);

	tx = 32; ty = 32; tz = 0;
	id = TOID(tx, ty, tz, texSize); 
	inputVoronoi[id] = ENCODE(tx, ty, tz); 
	inputPoints[nCount++] = ENCODE(tx, ty, tz);

	tx = 32; ty = 0; tz = 32;
	id = TOID(tx, ty, tz, texSize); 
	inputVoronoi[id] = ENCODE(tx, ty, tz); 
	inputPoints[nCount++] = ENCODE(tx, ty, tz);

	printf("nVertices: %d\n", nPoints);
	printf("nCount: %d\n", nCount);
}

void demo2(int texSize, int nPoints) {
	//int nVertices   = 2*2*2 + 7;
	int tx, ty, tz, id; 
	int nCount = 0;
	int offset = 10; 
	int step = 40;
	int size = 64;
	int len = (size - 2*offset) / step + 1;
	for(int i=0; i<len; i++) {
		tx = i*step + offset;
		for(int j=0; j<len; j++) {
			ty = j*step + offset;
			for(int k=0; k<len; k++) {
				tz = k*step + offset;
				id = TOID(tx, ty, tz, texSize); 
				inputVoronoi[id] = ENCODE(tx, ty, tz); 
				inputPoints[nCount++] = ENCODE(tx, ty, tz);
			}
		}
	}

	tx = 32; ty = 32; tz = 63;
	id = TOID(tx, ty, tz, texSize); 
	inputVoronoi[id] = ENCODE(tx, ty, tz); 
	inputPoints[nCount++] = ENCODE(tx, ty, tz);

	tx = 0; ty = 63; tz = 63;
	id = TOID(tx, ty, tz, texSize); 
	inputVoronoi[id] = ENCODE(tx, ty, tz); 
	inputPoints[nCount++] = ENCODE(tx, ty, tz);

	tx = 63; ty = 63; tz = 63;
	id = TOID(tx, ty, tz, texSize); 
	inputVoronoi[id] = ENCODE(tx, ty, tz); 
	inputPoints[nCount++] = ENCODE(tx, ty, tz);

	tx = 0; ty = 32; tz = 32;
	id = TOID(tx, ty, tz, texSize); 
	inputVoronoi[id] = ENCODE(tx, ty, tz); 
	inputPoints[nCount++] = ENCODE(tx, ty, tz);

	/*tx = 32; ty = 63; tz = 32;
	id = TOID(tx, ty, tz, texSize); 
	inputVoronoi[id] = ENCODE(tx, ty, tz); 
	inputPoints[nCount++] = ENCODE(tx, ty, tz);*/

	tx = 63; ty = 32; tz = 32;
	id = TOID(tx, ty, tz, texSize); 
	inputVoronoi[id] = ENCODE(tx, ty, tz); 
	inputPoints[nCount++] = ENCODE(tx, ty, tz);

	tx = 32; ty = 32; tz = 0;
	id = TOID(tx, ty, tz, texSize); 
	inputVoronoi[id] = ENCODE(tx, ty, tz); 
	inputPoints[nCount++] = ENCODE(tx, ty, tz);

	tx = 32; ty = 0; tz = 32;
	id = TOID(tx, ty, tz, texSize); 
	inputVoronoi[id] = ENCODE(tx, ty, tz); 
	inputPoints[nCount++] = ENCODE(tx, ty, tz);

	printf("nVertices: %d\n", nPoints);
	printf("nCount: %d\n", nCount);
}

void sphereDemo(int texSize, int nPoint) {
	int tx, ty, tz, id; 
	int nCount = 0;

	int cx = texSize/2, cy = texSize/2, cz = texSize/2;
	int offset = 5;
	int n = 5;

	tx = cx; ty = cy; tz = cz;
	id = TOID(tx, ty, tz, texSize); 
	inputVoronoi[id] = ENCODE(tx, ty, tz); 
	inputPoints[nCount++] = ENCODE(tx, ty, tz);
	double angle = 20.0;

	//6*n + 1
	for(int i=1; i<=1; i++) {
		tx = cx-10*i * sin(angle); ty = cy; tz = cz;
		id = TOID(tx, ty, tz, texSize); 
		inputVoronoi[id] = ENCODE(tx, ty, tz); 
		inputPoints[nCount++] = ENCODE(tx, ty, tz);

		tx = cx+10*i * sin(angle); ty = cy; tz = cz;
		id = TOID(tx, ty, tz, texSize); 
		inputVoronoi[id] = ENCODE(tx, ty, tz); 
		inputPoints[nCount++] = ENCODE(tx, ty, tz);

		tx = cx; ty = cy-10*i * sin(angle); tz = cz;
		id = TOID(tx, ty, tz, texSize); 
		inputVoronoi[id] = ENCODE(tx, ty, tz); 
		inputPoints[nCount++] = ENCODE(tx, ty, tz);

		tx = cx; ty = cy+10*i * sin(angle); tz = cz;
		id = TOID(tx, ty, tz, texSize); 
		inputVoronoi[id] = ENCODE(tx, ty, tz); 
		inputPoints[nCount++] = ENCODE(tx, ty, tz);

		tx = cx; ty = cy; tz = cz-10*i * sin(angle);
		id = TOID(tx, ty, tz, texSize); 
		inputVoronoi[id] = ENCODE(tx, ty, tz); 
		inputPoints[nCount++] = ENCODE(tx, ty, tz);

		tx = cx; ty = cy; tz = cz+10*i * sin(angle);
		id = TOID(tx, ty, tz, texSize); 
		inputVoronoi[id] = ENCODE(tx, ty, tz); 
		inputPoints[nCount++] = ENCODE(tx, ty, tz);
	} 

}

void demoTest(int texSize, int nPoints) {
	//int nVertices   = 2*2*2 + 9;
	int tx, ty, tz, id; 
	int nCount = 0;
	int offset = 10; 
	int step = 40;
	int len = (fboSize - 2*offset) / step + 1;
	for(int i=0; i<len; i++) {
		tx = i*step + offset;
		for(int j=0; j<len; j++) {
			ty = j*step + offset;
			for(int k=0; k<len; k++) {
				tz = k*step + offset;
				id = TOID(tx, ty, tz, texSize); 
				inputVoronoi[id] = ENCODE(tx, ty, tz); 
				inputPoints[nCount++] = ENCODE(tx, ty, tz);
			}
		}
	}

	tx = 32; ty = 32; tz = 63;
	id = TOID(tx, ty, tz, texSize); 
	inputVoronoi[id] = ENCODE(tx, ty, tz); 
	inputPoints[nCount++] = ENCODE(tx, ty, tz);

	tx = 0; ty = 63; tz = 40;
	id = TOID(tx, ty, tz, texSize); 
	inputVoronoi[id] = ENCODE(tx, ty, tz); 
	inputPoints[nCount++] = ENCODE(tx, ty, tz);

	tx = 63; ty = 63; tz = 63;
	id = TOID(tx, ty, tz, texSize); 
	inputVoronoi[id] = ENCODE(tx, ty, tz); 
	inputPoints[nCount++] = ENCODE(tx, ty, tz);

	tx = 0; ty = 32; tz = 32;
	id = TOID(tx, ty, tz, texSize); 
	inputVoronoi[id] = ENCODE(tx, ty, tz); 
	inputPoints[nCount++] = ENCODE(tx, ty, tz);

	tx = 32; ty = 63; tz = 32;
	id = TOID(tx, ty, tz, texSize); 
	inputVoronoi[id] = ENCODE(tx, ty, tz); 
	inputPoints[nCount++] = ENCODE(tx, ty, tz);

	tx = 63; ty = 32; tz = 32;
	id = TOID(tx, ty, tz, texSize); 
	inputVoronoi[id] = ENCODE(tx, ty, tz); 
	inputPoints[nCount++] = ENCODE(tx, ty, tz);

	tx = 32; ty = 32; tz = 0;
	id = TOID(tx, ty, tz, texSize); 
	inputVoronoi[id] = ENCODE(tx, ty, tz); 
	inputPoints[nCount++] = ENCODE(tx, ty, tz);

	tx = 32; ty = 0; tz = 32;
	id = TOID(tx, ty, tz, texSize); 
	inputVoronoi[id] = ENCODE(tx, ty, tz); 
	inputPoints[nCount++] = ENCODE(tx, ty, tz);

	tx = 15; ty = 15; tz = 50;
	id = TOID(tx, ty, tz, texSize); 
	inputVoronoi[id] = ENCODE(tx, ty, tz); 
	inputPoints[nCount++] = ENCODE(tx, ty, tz);


	printf("nVertices: %d\n", nPoints);
	printf("nCount: %d\n", nCount);
}




void cubeDemo(int texSize, int nPoints) {
	int len = 67;
	int step = 2;
	int nCount = 0;
	int tx, ty, tz, id;
	for(int i=0; i<len; i++) {
		tx = i*step;
		for(int j=0; j<len; j++) {
			ty = j*step;
			for(int k=0; k<len; k++) {
				tz = k*step;
				id = TOID(tx, ty, tz, texSize); 
				inputVoronoi[id] = ENCODE(tx, ty, tz); 
				inputPoints[nCount++] = ENCODE(tx, ty, tz);
			}
		}
	}
	printf("nCount: %d\n", nCount);
}

// Generate input points
void generateRandomPoints(int texSize, int nPoints)
{	
    for (int i = 0; i < texSize * texSize * texSize; i++)
        inputVoronoi[i] = MARKER; 
	
	//randomPoints(texSize, nPoints);
	demoFun(texSize, nPoints);
	//demo2(texSize, nPoints);

	//sphereDemo(texSize, nPoints);
}

/**********************************************************************************************
 * Deinitialization
 **********************************************************************************************/
void deinitialization()
{
    pba3DDeinitialization(); 

    free(inputPoints); 
    free(inputVoronoi); 
    free(outputVoronoi);
	if(tetrahedral != NULL)
		free(tetrahedral);

	for(int i=0; i<boundaryCount; i++) {
		if(boundaryMissing[i].tet != NULL)
			free(boundaryMissing[i].tet);
	}
	free(boundaryMissing);

}

/**********************************************************************************************
 * Initialization
 **********************************************************************************************/
bool initialization(char *fileName)
{
	printf("Read input file...\n"); 
	FILE *file = fopen(fileName, "r");
	if(file == NULL) {
		printf("can not access %s!!!\n", fileName);
		return false;
	}
	fscanf(file, "#%d %d %d\n", &fboSize, &nVertices, &display);

    pba3DInitialization(fboSize); 

    inputPoints     = (int *) malloc(nVertices * sizeof(int)); 
    inputVoronoi    = (int *) malloc(fboSize * fboSize * fboSize * sizeof(int)); 
    outputVoronoi   = (int *) malloc(fboSize * fboSize * fboSize * sizeof(int));

	//init the boundary missing tet
	boundaryMissing = (BoundaryMissing *)malloc(boundaryCount * sizeof(BoundaryMissing));
	for(int i=0; i<boundaryCount; i++) {
		boundaryMissing[i].tet = NULL;
		boundaryMissing[i].count = 0;
	}

	//read node data
	int id, tx, ty, tz, index;
	for (int i = 0; i < fboSize * fboSize * fboSize; i++)
		inputVoronoi[i] = MARKER;
	fscanf(file, "%d %d %d %d\n", &id, &tx, &ty, &tz);
	for(int i=0; i<nVertices; i++) {
		fscanf(file, "%d %d %d %d\n", &index, &tx, &ty, &tz);
		id = TOID(tx, ty, tz, fboSize); 
		inputVoronoi[id] = ENCODE(tx, ty, tz); 
		inputPoints[index] = ENCODE(tx, ty, tz);
	}

	fclose(file);
	return true;
	
}

/**********************************************************************************************
* Verify the output Voronoi Diagram
 **********************************************************************************************/
void compareResult(ErrorStatistics *e) 
{
    e->totalDistError = 0.0; 
    e->maxDistError = 0.0; 
    e->errorCount = 0; 

    int dx, dy, dz, nx, ny, nz; 
    double dist, myDist, correctDist, error;

    for (int i = 0; i < fboSize; i++)
        for (int j = 0; j < fboSize; j++) 
            for (int k = 0; k < fboSize; k++) {
                int id = TOID(i, j, k, fboSize); 
                DECODE(outputVoronoi[id], nx, ny, nz); 

                dx = nx - i; dy = ny - j; dz = nz - k; 
                correctDist = myDist = dx * dx + dy * dy + dz * dz; 

                for (int t = 0; t < nVertices; t++) {
                    DECODE(inputPoints[t], nx, ny, nz); 
                    dx = nx - i; dy = ny - j; dz = nz - k; 
                    dist = dx * dx + dy * dy + dz * dz; 

                    if (dist < correctDist)
                        correctDist = dist; 
                }

                if (correctDist != myDist) {
                    error = fabs(sqrt(myDist) - sqrt(correctDist)); 

                    e->errorCount++; 
                    e->totalDistError += error; 

                    if (error > e->maxDistError)
                        e->maxDistError = error; 
                }
            }
}

void printStatistics(ErrorStatistics *e)
{
    double avgDistError = e->totalDistError / e->errorCount; 

    if (e->errorCount == 0)
        avgDistError = 0.0; 

    printf("* Error count           : %i -> %.3f%\n", e->errorCount, 
        (double(e->errorCount) / nVertices) * 100.0);
    printf("* Max distance error    : %.5f\n", e->maxDistError);
    printf("* Average distance error: %.5f\n", avgDistError);
}

/**********************************************************************************************
 * Run the tests
 **********************************************************************************************/
void runTests()
{
    /*printf("Generate input...\n"); 
    generateRandomPoints(fboSize, nVertices); */

    //printf("Running PBA to compute 3D Voronoi Diagram...\n");
	printf("-----------------\n");
	printf("Texture: %dx%dx%d\n", fboSize, fboSize, fboSize);
	printf("Points: %d\n", nVertices);
	printf("-----------------\n");

	time_t start = clock();
	printf("========Start GPU DT 3D=======\n");
	pba3DVoronoiDiagram(inputVoronoi, outputVoronoi, phase1Band, phase2Band, phase3Band); 
	time_t end = clock();
	printf("========End GPU DT 3D run time: %fs==========\n", (double)(end - start) / CLOCKS_PER_SEC);
	//printf("inner tetraCount: %d\n", tetraCount);
	//printf("boundary tetreCount: ");
	int count = 0;
	for(int p=0; p<boundaryCount; p++) {		printf("%d ", boundaryMissing[p].count);
		count += boundaryMissing[p].count;
	}
	printf("\n");
	printf("total boundary tet count: %d\n", count);
	printf("total mesh tetrahedral : %d\n", count + tetraCount);

	//printf("Verifying the result...\n"); 
    //compareResult(&pba);

    //printStatistics(&pba); 
}

void generationTexture() {
	// Generate the color map
	GLubyte *colorMap = (GLubyte *) malloc(fboSize * fboSize * fboSize * 3 * sizeof(GLubyte)); 
	int tx, ty, tz;
	for (int i = 0; i < nVertices; i++) {
		GLubyte r = GLubyte(256 * random()); 
		GLubyte g = GLubyte(256 * random()); 
		GLubyte b = GLubyte(256 * random());

		DECODE(inputPoints[i], tx, ty, tz); 
		int id = TOID(tx, ty, tz, fboSize);

		colorMap[id * 3 + 0] = r; 
		colorMap[id * 3 + 1] = g;
		colorMap[id * 3 + 2] = b;
	}

	// Generate the voronoi diagram color 
	GLubyte *voronoi = (GLubyte *) malloc(fboSize * fboSize * fboSize * 3 * sizeof(GLubyte)); 

	for (int i = 0; i < fboSize; i++)
		for (int j = 0; j < fboSize; j++) 
			for (int k = 0; k < fboSize; k++) {
				int id = TOID(i, j, k, fboSize); 
				DECODE(outputVoronoi[id], tx, ty, tz);
				int cId = TOID(tx, ty, tz, fboSize);
				voronoi[id*3 + 0] = colorMap[cId*3 + 0];
				voronoi[id*3 + 1] = colorMap[cId*3 + 1];
				voronoi[id*3 + 2] = colorMap[cId*3 + 2];
			}

	/*for (int i = 0; i < 1; i++)
		for (int j = 0; j < fboSize; j++) 
			for (int k = 0; k < fboSize; k++) {
				int id = TOID(i, j, k, fboSize); 
				DECODE(outputVoronoi[id], tx, ty, tz);
				int cId = TOID(tx, ty, tz, fboSize);
				voronoi[id*3 + 0] = colorMap[cId*3 + 0];
				voronoi[id*3 + 1] = colorMap[cId*3 + 1];
				voronoi[id*3 + 2] = colorMap[cId*3 + 2];
			}*/


	/*for (int i = 0; i < fboSize; i++)
		for (int j = 0; j < fboSize; j++) 
			for (int k = fboSize-2; k <= fboSize-2; k++) {
				int id = TOID(i, j, k, fboSize); 
				int id2 = TOID(i, j, 0, fboSize);
				DECODE(nextBoundary[id2], tx, ty, tz);
				int cId = TOID(tx, ty, tz, fboSize);
				voronoi[id*3 + 0] = colorMap[cId*3 + 0];
				voronoi[id*3 + 1] = colorMap[cId*3 + 1];
				voronoi[id*3 + 2] = colorMap[cId*3 + 2];
			}*/

	//3D
	glPixelStorei(GL_UNPACK_ALIGNMENT, 1); 	texture = -1;	glGenTextures(1, &texture); 	glBindTexture(GL_TEXTURE_3D, texture); 	glTexImage3D(GL_TEXTURE_3D, 0, GL_RGB, fboSize, fboSize, 		fboSize, 0, GL_RGB, GL_UNSIGNED_BYTE, voronoi); 	glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE); 	glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE); 	glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_R, GL_CLAMP_TO_EDGE); 	glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MAG_FILTER, 		GL_LINEAR); 	glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MIN_FILTER, 		GL_LINEAR);

	free(voronoi);
	free(colorMap);
}

void glutDisplay() {
	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
	glLoadIdentity();
	glTranslatef(-0.5, -0.5f, -2.0f);
	if(rotateX < 0) rotateX += 360;
	if(rotateY < 0) rotateY += 360;
	glRotatef(rotateX, 1, 0, 0); 
	glRotatef(rotateY, 0, 1, 0); 
	glScalef(scale, scale, scale);

	glPushMatrix();
	glEnable(GL_DEPTH_TEST);
	//glDisable(GL_DEPTH_TEST);
		int t, tx, ty, tz;	if(faceColor == 1) {		glColor3f(0.2f, 0.8f, 0.2f);		for(int i=0; i<tetraCount; i++) {			for(int j=0; j<4; j++) {				glBegin(GL_TRIANGLES);				for(int k=j; k<j+3; k++) {					t = tetrahedral[i].vertice[k%4];					DECODE(t, tx, ty, tz);
					glVertex3f(tx, ty, tz);				}				glEnd();			}		}		glColor3f(0.2f, 0.8f, 0.2f);		for(int p=0; p<boundaryCount; p++) {			for(int i=0; i<boundaryMissing[p].count; i++) {				for(int j=0; j<4; j++) {					glBegin(GL_TRIANGLES);					for(int k=j; k<j+3; k++) {						t = boundaryMissing[p].tet[i].vertice[k%4];						DECODE(t, tx, ty, tz);
						glVertex3f(tx, ty, tz);					}					glEnd();				}			}		}	}	glLineWidth(3.0);	glColor3f(1.0f, 0.0f, 0.0f);	for(int p=0; p<boundaryCount; p++) {		for(int i=0; i<boundaryMissing[p].count; i++) {			for(int j=0; j<4; j++) {				glBegin(GL_LINE_LOOP);				for(int k=j; k<j+3; k++) {					t = boundaryMissing[p].tet[i].vertice[k%4];					DECODE(t, tx, ty, tz);
					glVertex3f(tx, ty, tz);				}				glEnd();			}		}	}		glColor3f(1.0f, 1.0f, 1.0f);	for(int i=0; i<tetraCount; i++) {		for(int j=0; j<4; j++) {			glBegin(GL_LINE_LOOP);			for(int k=j; k<j+3; k++) {				t = tetrahedral[i].vertice[k%4];				DECODE(t, tx, ty, tz);
				glVertex3f(tx, ty, tz);			}			glEnd();		}	}		/*glColor3f(0.2f, 0.8f, 0.2f);	glBegin(GL_TRIANGLES);		glVertex3f(0, 0, 0);		glVertex3f(0, 20, 0);		glVertex3f(20, 0, 0);			glEnd();	glColor3f(1.0f, 1.0f, 1.0f);	glBegin(GL_LINES);	glVertex3f(-10, 10, -10);	glVertex3f(30, 10, -10);	glEnd();*/	/*glColor3f(1.0f, 0.0f, 0.0f);	glBegin(GL_POINTS);		glVertex3f(0, 63, 63);	glEnd();*/	glColor3f(1.0f, 0, 0);	glBegin(GL_POINTS);		glVertex3f(0, 63, 63);	glEnd();	glColor3f(1.0f, 0, 0);	glBegin(GL_POINTS);	glVertex3f(63, 63, 63);	glEnd();		//glEnable(GL_DEPTH_TEST);	glPopMatrix();	glEndList();
	glutSwapBuffers();
}

void voronoiGlutDisplay() {
	glClear(GL_COLOR_BUFFER_BIT);
	glLoadIdentity();
	glTranslatef(0.0, 0.0f, -2.0f);
	if(rotateX < 0) rotateX += 360;
	if(rotateY < 0) rotateY += 360;
	glRotatef(rotateX, 1, 0, 0); 
	glRotatef(rotateY, 0, 1, 0); 
	glScalef(scale, scale, scale);

	glPushMatrix();
	//glDisable(GL_DEPTH_TEST);


	glEnable(GL_TEXTURE_3D);
	glBindTexture(GL_TEXTURE_3D,texture);	glBegin(GL_QUADS); 	for(int i = 0; i < 1000; i ++)	{		glTexCoord3f(0.0, 0.0, (GLfloat)i * 0.001); glVertex3f(-0.5, -0.5, (GLfloat)i * 0.001); 		glTexCoord3f(0.0, 1.0, (GLfloat)i * 0.001); glVertex3f(-0.5, 0.5, (GLfloat)i * 0.001); 		glTexCoord3f(1.0, 1.0, (GLfloat)i * 0.001); glVertex3f(0.5, 0.5, (GLfloat)i * 0.001); 		glTexCoord3f(1.0, 0.0, (GLfloat)i * 0.001); glVertex3f(0.5, -0.5, (GLfloat)i * 0.001); 	}	glEnd();		/*glColor3f(1.0f, 1.0f, 1.0f);	glBegin(GL_LINES);		glVertex3f(0, 0, 0);		glVertex3f(100, 0, 0);	glEnd();*/		//glEnable(GL_DEPTH_TEST);	glPopMatrix();	glEndList();
	glutSwapBuffers();
}


void resize(GLsizei w, GLsizei h) {
	glViewport(0,0,w,h);
	glMatrixMode(GL_PROJECTION);
	glLoadIdentity();
	gluPerspective(90.0f, (GLfloat)w/h, 0.01f, 100.0f);
	glMatrixMode(GL_MODELVIEW);
	glLoadIdentity();
}

void glut_mouse(int button, int state, int x, int y) 
{
	if (state == GLUT_UP)
		switch (button)
	{
		case GLUT_LEFT_BUTTON:
			isLeftMouseActive = false;
			break;
		case GLUT_RIGHT_BUTTON:
			isRightMouseActive = false; 
			break; 
	}

	if (state == GLUT_DOWN)
	{
		oldMouseX = x;
		oldMouseY = y;

		switch (button)
		{
		case GLUT_LEFT_BUTTON:
			isLeftMouseActive = true; 
			break;
		case GLUT_RIGHT_BUTTON:
			isRightMouseActive = true;
			break;
		}
	}		
}

void glut_mouseMotion(int x, int y) 
{
	if (isLeftMouseActive) {
		rotateY += 5.0 * double(x - oldMouseX) / scale / fboSize; 
		rotateX += 5.0 * double(y - oldMouseY) / scale / fboSize; 
		glutPostRedisplay(); 
	}
	else if (isRightMouseActive) {
		scale -= (y - oldMouseY) * scale / 400.0;
		glutPostRedisplay(); 
	} 

	oldMouseX = x; oldMouseY = y; 
}

int main(int argc, char **argv)
{
	char *fileName = "example.node";
	if (argc == 2) {
		fileName = argv[1];
	}
	glutInitWindowPosition(0, 0); 
	glutInitWindowSize(600, 600); 
	glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE | GLUT_DEPTH | GLUT_ALPHA); 
	glutInit(&argc, argv); 
	glutCreateWindow("3D PBA");
	glewInit();

	if(! initialization(fileName))
		return 0;
     
    runTests(); 

	/*generationTexture(); 
	scale = 1.0f;
	glutDisplayFunc(voronoiGlutDisplay); */

	if(display) {
		glutDisplayFunc(glutDisplay); 
		glutReshapeFunc(resize);
		glutMouseFunc(glut_mouse); 
		glutMotionFunc(glut_mouseMotion); 
		glutMainLoop();
	}

	deinitialization(); 

    getchar(); 
		
	return 0;
}